Heatmap

Preparation

Paths and libraries setting

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/paths.R")
source("../../../source/functions.R")

# Florentin CONSTANCIAS's script
source("https://raw.githubusercontent.com/fconstancias/DivComAnalyses/master/R/phyloseq_heatmap.R")

Load phyloseq objects before and after decontam

setwd(path_rdata)
ps <- readRDS("MED_phyloseq.rds")
ps.filter <- readRDS("1D_MED_phyloseq_decontam.rds")

Change factors format for organ and species

change_organ(ps.filter)
change_species(ps.filter)

Heatmap

All the MED nodes

physeq <- ps.filter
colnames(tax_table(physeq))[7] <- "Strain"

levels(physeq@sam_data$Species) <- c("AA", "CP", "CQ")
# levels(physeq@sam_data$Species)

levels(physeq@sam_data$Strain) <- c("B", "CE", "G", "L", "W-")
# levels(physeq@sam_data$Strain)

physeq %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Genus",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Loading required package: ampvis2
## Warning: There are only 13 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free")+
  scale_fill_viridis_c(breaks = c(0,  0.01, 0.1, 1, 10, 100), 
                       labels = c(0,  0.01, 0.1, 1, 10, 100), 
                      trans = scales::pseudo_log_trans(sigma=0.001), # Scaling factor for the linear part of pseudo-log transformation
                       na.value = 'transparent') -> p1
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p1[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p1

Wolbachia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Wolbachia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 20 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p2
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p2[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p2

Asaia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Asaia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 24 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p3
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p3[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p3

Elizabethkingia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Elizabethkingia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 5 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p4
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p4[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p4

Legionella nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Legionella") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 3 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p5
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p5[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p5

Chryseobacterium nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Chryseobacterium") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 4 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p6
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p6[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p6

Erwinia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Erwinia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 3 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p7
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p7[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p7

Serratia nodes

physeq %>%
 transform_sample_counts(function(x) x/sum(x) *100) %>% 
 subset_taxa(Genus == "Serratia") %>%
 phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Strain", "Organ", "Field", "Species" ),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
## Warning: There are only 2 taxa, showing all
p + facet_grid( ~ Species + Strain + Organ + Field , scales = "free", space = "free") + 
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p8
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p8[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p8

Check Blanks

# import phyloseq object before decontam
setwd(path_rdata)
ps <- readRDS("MED_phyloseq.rds")
colnames(tax_table(ps))[7] <- "Strain"

# 70 nodes (all)
ps %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>% 
   phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Control"),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p
p + facet_grid( ~ Control, scales = "free", space = "free") +
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p9
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p9[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p9

# Pseudomonas nodes (3)
ps %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>% 
  subset_taxa(Genus=="Pseudomonas") %>%
   phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Control"),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p10
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p10[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p10

# Enhydrobacter nodes (3)
ps %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>% 
  subset_taxa(Genus=="Enhydrobacter") %>%
   phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Control"),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
    theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p11
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p11[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p11

# Rahnella1 nodes (3)
ps %>%
  transform_sample_counts(function(x) x/sum(x) *100) %>% 
  subset_taxa(Genus=="Rahnella1") %>%
   phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "SampleID",
                          facet_by = c("Control"),
                          tax_aggregate = "Species",
                          tax_add = NULL,
                          ntax  = 100) -> p_bis
p_bis + facet_grid( ~ Control, scales = "free", space = "free") +
  theme(axis.text.y = element_text(face="italic", angle = 0,  size = 12))+
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 100), 
                       labels = c(0,  0.01, 1, 10, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.001),
                       na.value = 'transparent') -> p12
## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.
p12[["labels"]][["fill"]] <- "Read abundance \n(pseudo log trans sigma=0.001)"
p12

Save heatmaps

setwd(path_plot)

tiff("1Ed_MED_heatmap.tiff", units="in", width=28, height=16, res=300)
p1
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_wolbachia.tiff", units="in", width=22, height=10, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_asaia.tiff", units="in", width=22, height=10, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_elizabethkingia.tiff", units="in", width=22, height=10, res=300)
p4
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_legionella.tiff", units="in", width=22, height=10, res=300)
p5
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_chryseobacterium.tiff", units="in", width=22, height=10, res=300)
p6
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_erwinia.tiff", units="in", width=22, height=10, res=300)
p7
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_serratia.tiff", units="in", width=22, height=10, res=300)
p8
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_blanks_all.tiff", units="in", width=22, height=10, res=300)
p9
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_blanks_pseudomonas.tiff", units="in", width=22, height=10, res=300)
p10
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_blanks_enhydrobacter.tiff", units="in", width=22, height=10, res=300)
p11
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ed_MED_heatmap_blanks_rahnella1.tiff", units="in", width=22, height=10, res=300)
p12
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap.png", units="in", width=22, height=10, res=300)
p1
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_wolbachia.png", units="in", width=22, height=10, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_asaia.png", units="in", width=22, height=10, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_elizabethkingia.png", units="in", width=22, height=10, res=300)
p4
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_legionella.png", units="in", width=22, height=10, res=300)
p5
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_chryseobacterium.png", units="in", width=22, height=10, res=300)
p6
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_erwinia.png", units="in", width=22, height=10, res=300)
p7
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_serratia.png", units="in", width=22, height=10, res=300)
p8
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_blanks_all.png", units="in", width=22, height=10, res=300)
p9
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_blanks_pseudomonas.png", units="in", width=22, height=10, res=300)
p10
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_blanks_enhydrobacter.png", units="in", width=22, height=10, res=300)
p11
dev.off()
## quartz_off_screen 
##                 2
png("1Ed_MED_heatmap_blanks_rahnella1.png", units="in", width=22, height=10, res=300)
p12
dev.off()
## quartz_off_screen 
##                 2